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Abstract. Previous work in the literature has studied gravitational radiation in black-hole 
collisions at the speed of light. In particular, it had been proved that the perturbative field 
equations may all be reduced to equations in only two independent variables, by virtue of 
a conformal symmetry at each order in perturbation theory. The Green function for the 
perturbative field equations is here analyzed by studying the corresponding second-order 
hyperbolic operator with variable coefficients, instead of using the reduction method from 
the retarded flat-space Green function in four dimensions. After reduction to canonical 
form of this hyperbolic operator, the integral representation of the solution in terms of the 
Riemann function is obtained. The Riemann function solves a characteristic initial-value 
problem for which analytic formulae leading to the numerical solution are derived. 
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1. Introduction 



The construction of suitable inverses of differential operators lies still at the very heart of 
many profound properties in classical and quantum field theory. For example, the theory 
of small disturbances in local field theory can only be built if suitable invertible operators 
are considered [1] . In a path-integral formulation, these correspond to the gauge-field and 
ghost operators, respectively [2, 3]. Moreover, the Peierls bracket on the space of physical 
observables, which is a Poisson bracket preserving the invariance under the full infinite- 
dimensional symmetry group of the theory, is obtained from the advanced and retarded 
Green functions of the theory via the supercommutator function [1-4] , and leads possibly 
to a deeper approach to quantization. Last, but not least, a perturbation approach to 
classical general relativity relies heavily on a careful construction of Green functions of 
operators of hyperbolic [5-8] and elliptic [9, 10] type. In particular, following [5-8], we 
shall be concerned with the axisymmetric collision of two black holes travelling at the speed 
of light, each described in the centre-of-mass frame before the collision by an impulsive 
plane-fronted shock wave with energy [i. One then passes to a new frame to which a large 
Lorentz boost is applied. There the energy v = \ie a of the incoming shock 1 obeys v » A, 

where A = [ie~ a is the energy of the incoming shock 2 and e Q = yj j^jj ((3 being the usual 

relativistic parameter). In the boosted frame, to the future of the strong shock 1, the 
metric can be expanded in the form [6, 8] 



L i=l v 7 J 

where r\ a b is the standard notation for the Minkowski metric. The task of solving the 
Einstein field equations becomes then a problem in singular perturbation theory, having 



spectively in ^, once that characteristic initial data are given just to the future of the 
strong shock 1. The perturbation series (1.1) is physically relevant because, on boosting 
back to the centre-of-mass frame, it is found to give an accurate description of space-time 
geometry where gravitational radiation propagates at small angles away from the forward 




(1.1) 



to find , ... by solving the linearized field equations at first, second, ... order re- 
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symmetry axis 9 = 0. The news function cq, which describes gravitational radiation arriv- 
ing at future null infinity in the centre-of-mass frame, is expected to have the convergent 
series expansion [6, 8] 



with f a suitable retarded time coordinate, and \i the energy of each incoming black hole in 
the centre-of-mass frame. In [6, 8] a very useful analytic expression of aiirjn) was derived, 
exploiting the property that perturbative field equations may all be reduced to equations 
in only two independent variables, by virtue of a remarkable conformal symmetry at each 
order in perturbation theory. The Green function for perturbative field equations was then 
found by reduction from the retarded flat-space Green function in four dimensions. 

However, a direct approach to the evaluation of Green functions appears both desirable 
and helpful in general, and it has been our aim to pursue such a line of investigation. 
For this purpose, reduction to two dimensions with the associated hyperbolic operator 
is studied again in section 2. Section 3 performs reduction to canonical form with the 
associated Riemann function. Equations for the Goursat problem obeyed by the Riemann 
function are derived in section 4, while the corresponding numerical algorithm is discussed 
in section 5. 

2. Reduction to two dimensions and the associated operator 

As is well known from the work in [6] and [8], the field equations for the first-order cor- 
rection i n the expansion (1.1) are particular cases of the general system given by the 
flat-space wave equation (here u = -h=(z + t),v = -ys(z — t)) 



oo 




(1.2) 



n=0 




(2.1) 



supplemented by the boundary condition 



^{u = 0) = e im< V- n .f[81ogM - V2v] 



(2.2a) 
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/Or) =0 Vx < 0. (2.26) 
Moreover, ip should be of the form e ^m ^p _n x(q , , r) for u > 0, where 

g = up~ 2 (2.3) 

r = 81og(^p) - V2v. (2.4) 

For the homogeneous wave equation (2.1) there is no advantage in eliminating p and <p 
from the differential equation. However, the higher-order metric perturbations turn out to 
obey inhomogeneous flat-space wave equations of the form 

□ V = 5 (2.5) 

where £ is a source term equal to e im ^ p~( n + 2 ) H(q, r). This leads to the following equation 

for x = e" im0 p n ip: 

£>m,nX(q,r) = H(q,r) (2.6) 

where C m ,n is an hyperbolic operator in the independent variables q and r, and takes the 
form [6, 8] 

02 ^2 ^2 

£ m , n = -(2^ + 32^)— — + 4q 2 — + 64 ° 



dqdr dq 2 dr 2 

+ 4(n + l)q— -16n— +n 2 -m 2 . (2.7) 
oq or 

The proof of hyperbolicity of £ m , n , with the associated normal hyperbolic form, can be 
found in section 3 of [6], and in [8]. The advantage of studying Eq. (2.6) is twofold: to 
evaluate the solution at some space-time point one has simply to integrate the product of 
H and the Green function G m ^ n of C m ,n- 



X(q, r) = J G m , n (q, r; q , r )H(q , ro)dq Q dr (2.8) 
and the resulting numerical calculation of the solution is now feasible [7, 8]. 
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Since we are interested in a direct approach to the evaluation of the Green function 
G m , n in the (q, r) coordinates, we begin by noticing that, in all derivatives which are not 
mixed, the operator C m ,n can be made a constant coefficient operator upon setting 



which implies 



a = log(g) (2.9) 



and hence 



2 d 2 d 2 d_ 

^ dq 2 da 2 da 



q2 q2 q2 

C m , n = -(2V2e-« + 32)-— + 4— + 64- 



This operator is further simplified upon defining the variable 



4 

which implies that 



1 \ d 2 ( d 2 d 2 



so that 

d 2 d 2 d 2 



(2.11) 



dadr da 2 dr 2 

+ 4n- 16ra— + n 2 -m 2 . (2.12) 

da or 



R=- A (2.13) 



y/2 J dadR \da 2 dR 2 

+ 4n(- ^)+™ 2 -™ 2 - (2-14) 



<9ct <9.R/ 
This suggests defining yet new variables 

X = a + R (2.15) 
Y = a-R (2.16) 



dadtf ^ 2 dY 2 



(2.17) 
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*L + *- = 2(*-+*-\ (218) 
da 2 8R 2 \dx 2 dY 2 ) K ' 

d d d , 



and hence C m ,n reads eventually 

T m , n = 16^j + Sn A + n 2 - m 2 

with Green function satisfying the equation 

T m , n G m , n (e 2 ^ , 2(X - Y); e^, 2(X - Y )) 

= \& (e 2 ^ - e™ ) - y) - (X - Y )). (2.21) 

The operator T m n is the sum of an elliptic operator in the Y variable and a two-dimensional 
wave operator 'weighted' with the exponential e~( x+Y ^ 2 , which is the main source of 
technical complications in these variables. 

3. Reduction to canonical form and the Riemann function 

It is therefore more convenient, in our general analysis, to reduce first Eq. (2.6) to canonical 
form, and then find an integral representation of the solution. Reduction to canonical 
form means that new coordinates x = x(q, r) and y = y(q, r) are introduced such that the 
coefficients of and Jp- vanish. As is shown in [6, 8], this is achieved if 



3x dy 

dr dr ^ ^ 



dx l + 8qV2 + A/1 + 16QV2 



dq 2V2q 2 



dy = l + 8qV2- y/l + 16qy/2 

dq 2V2q 2 ' 1 ' ° J 
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The resulting formulae are considerably simplified if one defines 



t = \Jl + 16qV2 = t(x,y). (3.4) 
The dependence of t on x and y is obtained implicitly by solving the system [6, 8] 

I = r + log(*zi)-^l ly -4 (3.5) 
!, = r + log(^) + ^-4. (3.6) 

This leads to the equation 

i (t - 1) 2t (x-y) , n „ . 

which can be cast in the form 



(t — 1) 2t A 
- -<=> (i 

+ 



e(^)= e V. (3.76) 



This suggests defining 



so that one first has to solve the transcendental equation 

we 2 » = e s (3.9) 
to obtain w = w(x — y), from which one gets 

t=|^| = ^-y). (3.10) 

On denoting by (7(10) the left-hand side of Eq. (3.9), one finds that, in the plane (w, g(w)), 
the right-hand side of Eq. (3.9) is a line parallel to the w-axis, which intersects g(w) at 
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no more than one point for each value of x — y. For example, when w = 1, g{w) intersects 
the line taking the constant value 1, for which x — y = 0. The function 



g : w — > g(w) = we 
is asymmetric and has the limiting behaviour described by 

lim g(w) = — oo lim g(w) = (3-H) 

lim g(w) = lim g(w) = oo. (3-12) 

Thus, in the lower half-plane, g has an horizontal asymptote given by the w-axis, and 
a vertical asymptote given by the line w = 0, while it has no asymptotes in the upper 
half-plane, since 

lim = oo 

in addition to (3.12). The first derivative of g reads 

g>) = <^±)!e^. (3.13) 

One therefore has g'(w) > for all w > 0, and g'(w) < for all w E (— oo, 0) — {—1}, and 
g is monotonically decreasing for negative w and monotonically increasing for positive w. 
The point w = —1, at which g'(w) vanishes, is neither a maximum nor a minimum point, 
because 

These formulae imply that <7"(— 1) =0 but g"'(— 1) = —1^0, and hence = — 1 yields a 
flex of <7(iy). 
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In the (x, y) variables, the operator C m ,n reads therefore 



d 2 d d 

£m,n = f( x 'y)^fry + 9{x,y)— + h(x,y)— + n 2 - m 2 



where, exploiting the formulae 



dx 64\/2 



dq (t-1) 2 

dy 64\/2 
dq~ ~ (£ + l) 2 



one finds 



,(,,,) = - (2 V2 + 3 25 )(| + |) + ^| + 12 8 

2t 2 (t 2 + 1) 



= 256 



dx 



g(x, y) = 4(n + 1)<7 ^ 16n = 16 

/i(>, y) = 4(n + - 16n = 16 



1 + 



2(n + l) 
(t-1) 



1 2(n + l) 



The resulting canonical form of Eq. (2.6) is 

L[X] = {^y + o( *' v) Tx + b{x > y) Yy + C(X ' V) ) X(X ' V) 
= H(x,y) 



where 



a(x,y) = 
b(x,y) = 
c(x,y) = 



g(x,y) = l (l-t)(t + l) 2 (2n + l + t) 

f(x,y) 16 (t 4 + 4t 2 -l) 

h(x,y) 1 (t + l)(t-l) 2 (2n + l-t) 
f(x,y) "16 (t 4 + 4t 2 -l) 

n 2 - m 2 _ (m 2 - n 2 ) (t - l) 2 (t + l) 2 



f(x,y) 



256 (t 4 + 4t 2 - 1) 
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H(x,y) = 



H(x,y) H(x,y)(t-l) 2 (t + lf 



(3.24) 



f(x,y) 256 (t 4 + 4t 2 -l) ' 

Note that a(-t) = b(t),b(-t) = a(t), c(-t) = c(t),H(-t) = H(t). 

For an hyperbolic equation in the form (3.20), we can use the Riemann integral rep- 
resentation of the solution. For this purpose, recall from [11] that, on denoting by the 
adjoint of the operator L in (3.20), which acts according to 

L] \x] = Xxy ~ (ax)x ~ (bx)y + cx (3-25) 

one has to find a 'function' R(x,y;£,rj) (actually a kernel) subject to the following condi- 
tions ((£,77) being the coordinates of a point P such that characteristics through it intersect 
a curve C at points A and P, AP being a segment with constant y, and BP being a segment 
with constant x): 



(i) As a function of x and y, R satisfies the adjoint equation 

(ii) R x = bR on AP, i.e. 

Rx(x, y; f , 77) = b(x, rj)R{x, y; f , 77) on y = 77 
and P y = aP on PP, i.e. 

#v(a, V\ ^ V) = y)R(x, y; f , 77) on x = £ 
(hi) P equals 1 at P, i.e. 

P(£,77;£,?7) = l. 

It is then possible to express the solution of Eq. (3.20) in the form 

x{P) = \[x{A)R{A) + x{B)R{B)] 



+ 



^\., + (bR - X -R x ) \ 



AB 

7T\.</ + [aR-^R, ) \ 



dx 



dy 



+ 



/ / R(x,y;£,rj)H(x,y)dxdy 

J JQ, 



(3.26) 



(3.27) 



(3.28) 



(3.29) 



(3.30) 
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where Q is a domain with boundary. 

Note that Eqs. (3.27) and (3.28) are ordinary differential equations for the Riemann 
function R(x, y; £, rj) along the characteristics parallel to the coordinate axes. By virtue of 
(3.29), their integration yields 



which are the values of R along the characteristics through P. Equation (3.30) yields 
instead the solution of Eq. (3.20) for arbitrary initial values given along an arbitrary 
non-characteristic curve C, by means of a solution R of the adjoint equation (3.26) which 
depends on x, y and two parameters £, rj. Unlike the Riemann function R solves a 
characteristic initial-value problem. 

4. Goursat problem for the Riemann function 

By fully exploiting the reduction to canonical form of Eq. (2.6) we have considered novel 
features with respect to the analysis in [6, 8], because the Riemann formula (3.30) con- 
tains also the integral along the piece of curve C from A to £?, and the term ^[x(A)R(A) + 
x(B)R(B)]. This representation of the solution might be more appropriate for the numer- 
ical purposes considered in [7], but the task of finding the Riemann function R remains 
extremely difficult. One can however use approximate methods for solving Eq. (3.26). For 
this purpose, we first point out that, by virtue of Eq. (3.25), Eq. (3.26) is a canonical 
hyperbolic equation of the form 




(3.31) 




(3.32) 




(4.1) 



where 



A = - 



a 



(4.2) 



B 



= -b 



(4.3) 
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C = c - a x - by. (4.4) 

Thus, on denning 

U = R (4.5) 

V = R X + BR (4.6) 

the equation (4.1) for the Riemann function is equivalent to the hyperbolic canonical 
system [11] 

U x = fi(x,y)U + f 2 (x,y)V (4.7) 



where 



V y = gi {x,y)U + g 2 {x,y)V (4.8) 

h = -B = b (4.9) 

h = 1 (4-10) 

gi = AB - C + By = ab - c + a x (4.11) 

g 2 = -A = a. (4.12) 

For the system described by Eqs. (4.7) and (4.8) with boundary data (3.31) and (3.32) an 
existence and uniqueness theorem holds (see [11] for the Lipschitz conditions on boundary 
data), and we can therefore exploit the finite differences method to find approximate 
solutions for the Riemann function R(x,y;£,r]), and eventually x(-P) with the help of the 
integral representation (3.30). 

5. Concluding remarks 

The inverses of hyperbolic operators [12] and the Cauchy problem for hyperbolic equations 
with polynomial coefficients [13] have always been the object of intensive investigation in 
the mathematical literature. We have here considered the application of such issues to 
axisymmetric black hole collisions at the speed of light, relying on the work in [5-8]. We 
have pointed out that, for the inhomogeneous equations (2.6) occurring in the perturbative 
analysis, the task of inverting the operator (2.7) can be accomplished with the help of the 
Riemann integral representation (3.30), after solving Eq. (4.1) for the Riemann function. 
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One has then to solve a characteristic initial-value problem for a homogeneous hyperbolic 
equation in canonical form in two independent variables, for which we have developed 
formulae to be used for the numerical solution with the help of a finite differences scheme. 
For this purpose one studies the canonical system (cf (4.7) and (4.8)) 

U x = F(x,y,U,V) (5.1) 

V y = G(x,y,U,V) (5.2) 

in the rectangle 1Z = {x, y : x G [xq, xq + a], y G [yo, yo + 6]} with known values of U on the 
vertical side AD where x = xq, and known values of V on the horizontal side AB where 
y = yo. The segments AB and AD are then divided into m and n equal parts, respectively. 
On setting ^ = h and ^ = k, the original differential equations become equations relating 
values of U and V at three intersection points of the resulting lattice, i.e. 

u(p,., s+l) -u(p,. s) = F (8 3o) 

k 

It is now convenient to set U rs = U(P rs ), V rs = V(P rs ), so that these equations read 

U ri8+1 = U rs + hF(P rs , U rs , V rs ) (5.36) 

Vr+i,* = V rs + kG(P rs , U rs , V rs ). (5.46) 

Thus, if both U and V are known at P rs , one can evaluate U at P r)S +i and V at P r +i jS - 
The evaluation at subsequent intersection points of the lattice goes on along horizontal or 
vertical segments. In the former case, the resulting algorithm is 

s-l 



Urs =U r0 + hJ2 F(P rl , U ri , V ri ) (5.5) 
i=l 

V rs = V r -l,s + kG(P r -l,s, U r - lt8 , V r -l,s) (5.6) 
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while in the latter case one obtains the algorithm expressed by the equations 



r-1 



Vrs — Vqs + ^2 ^(-^S' Uisi Vis) 



(5.7) 



U rs — U r}S -l + hF(P rjS -i, U r>8 -1, V r}S -l). 



(5.8) 



Stability of such solutions is closely linked with the geometry of the associated character- 
istics, and the criteria to be fulfilled are studied in section 13.2 of [14] (stability depends 
crucially on whether or not | < 1). 

To sum up, one solves numerically Eq. (3.9) for w = w(x, y) = w(x — y), from which 
one gets t(x — y) with the help of (3.10), which is a fractional linear transformation. This 
yields a, 6, c and H as functions of (x,y) according to (3.21)-(3.24), and hence A, B and 
C in the equation for the Riemann function are obtained according to (4.2)-(4.4), where 
derivatives with respect to x and y are evaluated numerically. Eventually, the system given 
by (4.7) and (4.8) is solved according to the finite- differences scheme of the present section, 
with 



Once the Riemann function R = U is obtained with the desired accuracy, numerical 
evaluation of the integral (3.30) yields x(-P), and x((/,r) is obtained upon using Eqs. (3.5) 
and (3.6) for the characteristic coordinates. Our steps are conceptually desirable since 
they rely on well established techniques for the solution of hyperbolic equations in two 
independent variables [11, 14], and provide a viable alternative to the numerical analysis 
performed in [7], because all functions should be evaluated numerically. Our method is 
not obviously more powerful than the one used in [5-8] , but is well suited for a systematic 
and lengthy numerical analysis, while its analytic side provides an interesting alternative 
for the evaluation of Green functions both in black hole physics and in other problems 
where hyperbolic operators with variable coefficients might occur. This task remains very 
important because a strong production of gravitational radiation is mainly expected in the 



F = hU + f 2 V = hR + f 2 (R x + BR) 



(5.9) 



G = 9l U + g 2 V = 9l R + g 2 (R x + BR). 



(5.10) 
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extreme events studied in [5-8] and which motivated our paper. Any viable way of looking 
at mathematical and numerical aspects of the problem is therefore of physical interest for 
research planned in the years to come [15]. 
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